Matrix controlled channel diffusion of sodium in amorphous silica 
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To find the origin of the diffusion channels observed in 
sodium-silicate glasses, we have performed classical molecu- 
lar dynamics simulations of Na2 0-4Si02 during which the 
mass of the Si and O atoms has been multiplied by a tuning 
coefficient. We observe that the channels disappear and that 
the diffusive motion of the sodium atoms vanishes if this coef- 
ficient is larger than a threshold value. Above this threshold 
the vibrational states of the matrix are not compatible with 
those of the sodium ions. We interpret hence the decrease of 
the difi'usion by the absence of resonance conditions. 

PACS numbers: 61.43.Fs, 61.43.Bn, 66.30.Hs, 63.50.-Fx 



The mechanism of ionic transport in amorphous ma- 
terials [1,2] at the atomic scale is not completely eluci- 
dated so far and therefore it is the subject of numerous 
experimental [3,4] and numerical [5-9] studies. Among 
the most studied systems are the sodium-silicate glasses 
since they contain the essential ingredients, namely the 
amorphous matrix (silica) and the mobile ions (sodium), 
as a first step in the simulation of more complex glasses 
of higher practical interest. 

In previous studies [10-12], we have shown by means of 
classical molecular dynamics simulations of Na20-4Si02 
(NS4) , that the sodium atoms diffuse through a well con- 
nected network of pockets (which represents only a lim- 
ited fraction of the entire available space) that we have 
called "channels" to be coherent with the literature. The 
existence of the channels, which are not due to micro- 
segregation effects [6,7], has been confirmed by the exis- 
tence of a pre-peak in the partial Na-Na structure fac- 
tor at a wave- vector q = 0.95 A^^ [13]. This pre-peak 
has also been observed experimentally [14] and numer- 
ically [15] in another study. We have also shown that 
the location of the channels is strongly correlated to the 
positions of the non-bridging oxygens [13] and Horbach 
et at have shown that the sodium dynamics should be 
related to that of the underlying sihca network [15]. This 
suggests that the origin of the channels could be related 
to the dynamical properties of the matrix. To check this 
idea we present in this letter classical molecular dynam- 
ics simulations on a series of "toy" systems in which the 
atomic masses of both the oxygen and silicon atoms have 
been systematically changed after artificially multiplying 
their experimental values by a common factor varying 



from 0.5 to infinity, the usual NS4 system being recovered 
for /i = 1. 

By studying both the mean square displacement of the 
sodium atoms and the characteristics of the channels, we 
find a change in the sodium diffusion properties when 
the parameter fi is increased above a value of about 30. 
Above this threshold, the sodium diffusion decreases and 
the channels can no more be clearly defined. Guided by 
the concomitant change in the short time characteristics 
of the velocity autocorrelation function, we have calcu- 
lated the vibrational density of states (VDOS) for both 
the sodium atoms and the atoms of the matrix for various 
values of /i. We observe that the threshold corresponds 
to the value of fi above which the sharp VDOS of the 
sodium atoms starts to separate from the band of vibra- 
tional states of the matrix. Therefore we propose that the 
diffusive motion of the sodium atoms, i.e. their ability 
to escape from their local cage, is facilitated when their 
vibrational frequencies are compatible with the ones of 
the matrix. We argue that this mechanism is responsible 
for the channel diffusion of sodium in the silica matrix. 

In this letter we present classical molecular dynamics 
calculations of a system of 648 particles (86 sodium, 173 
silicon and 389 oxygen atoms) confined in a cubic box of 
edge length 20.88 A with periodic boundary conditions. 
The density is thus the experimental density of glassy 
NS4, i.e. 2.38 g.cm^'^ [16]. The interactions between 
the particles are given by a modified version of the so- 
called "BKS" potential [17,18] which is able to reproduce 
the structure as well as the dynamics of several sodium- 
silicate systems [10-13,19] (for more details see [11]). In 
the present study we have generated three independent 
samples (in the following, all the results are averaged 
over these samples in order to improve the statistics) 
at a temperature of ~ 1900 K for which the channels 
have been evidenced in our previous studies. For each 
sample we have performed ten different simulations in 
which the mass of the atoms of the matrix has been arti- 
ficially multiplied by a factor fj. of 0.5, 1, 2, 5, 10, 30, 10^, 
lO'^, 10^, and 10^. We have also carried out the limiting 
case fj, = oo hy performing simulations in which only the 
sodium atoms move, the atoms of the matrix being kept 
fixed (frozen matrix approximation). Each of these 33 
samples has first been relaxed at ^ 1900 K for 10® steps 
(1 step = 1.4 fs) and the following 10® steps were used 
to produce 1000 configurations recorded every 1000 time 
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steps. These configurations have then been used to an- 
alyze the trajectories of the sodium atoms during these 
1.4 ns. 

In Fig. 1 wc show B^{t), the mean square displace- 
ment (MSD) of the sodium atoms, for various values of 
/i. While these curves have the characteristic long-time 
shape of strong diffusion for small values of /i, they flat- 
ten out as j2 increases. In the limiting case 11 = 00 (frozen 
matrix) the curve becomes even completely flat, at least 
within the period of time covered by our simulations. In 
order to give a quantitative idea of how the ionic diffu- 
sion properties decrease with increasing /x, we have deter- 
mined a characteristic time tmsd necessary for B?{t) to 
be equal to 4 A^. This corresponds to an average travel 
distance of 2 A which is the minimum distance between 
two Na neighbors (as determined from the Na-Na radial 
pair distribution function [12]). In the inset of Fig. 1 we 
give the variation of tmsd with ji. While tmsd is almost 
independent of /i for small /i, it starts to increase for 
values of /x larger than a value of about 30. 

Since we have previously shown that the difl[usive 
motion of the sodium atoms takes place within chan- 
nels [10-13], it is worth studying how these channels are 
modified when the mass of the atoms of the matrix is 
changed. Therefore, for all our values of /i, we have 
determined the channels, in the same manner as pre- 
viously [10]: using a 3-dimensional mesh we determine 
the number of different sodium atoms that have visited 
each cube of the mesh during the total simulation time. 
Then we define ^ which is the minimal occupation num- 
ber such that the cubes visited more than ^ times rep- 
resent the upper 10 % of all the visited cubes (for more 
details see [10]). In previous studies we have shown that 
at T ^ 1900 K a cube needed to be visited by at least 
^ = 8 different sodium atoms during a 1.4 ns simulation, 
in order to be part of the channel structure. 

Here we stiidy how ^ changes with /i and the results 
(averaged over three samples) are reported in Fig. 2. 
While for fi < 30, ^ remains almost constant between 7 
and 8 within the error bars, it starts to decrease dramat- 
ically for fi > 30. Obviously, as explained in refs. [10 -13], 
one can no more speak of "channels" if ^ becomes as small 
as unity. Therefore one can interpret the results in Fig. 2 
by assuming that the channels disappear progressively as 
soon as n > 30. Then one can argue that the increase 
of the characteristic time for the diffusion observed in 
Fig. 1 is intimately correlated with the disappearance of 
the channels shown in Fig. 2. 

In order to know if the observed change in the long 
time diffusive behavior of the sodium atoms is accompa- 
nied by a change in their short time dynamical proper- 
ties, we have calculated the Na velocity autocorrelation 
function ^{t) = {v{to)v{to + 1)) / (v'^{to)) and studied its 
short time behavior. As seen in Fig. 3 the short time part 
of '&{t) is characterized by a first minimum at r^, typical 
of the time for a sodium atom to bounce back and forth 



against the internal boundaries of the cage in which it vi- 
brates (well before being eventually able to jump towards 
another cage). In the inset of Fig. 3 the intensity of the 
minimum 'i^m = 'i?(r^) has been plotted as a function of 
11. f^ni — 0.28 for /X < 10 while drops down to —0.5 for 
/i > 30. This increasing depth can be interpreted as an 
increasing stiffness of the internal cage boundaries when 
/X is increased. In a classical mechanics picture, due to 
this increased stiffness, the probability for a sodium atom 
to escape from its cage becomes weaker and this is con- 
sistent with the vanishing diffusion observed in Fig. 1. 

To go further in our microscopical analysis, and to 
understand the changes occurring around fi ^ 30, we 
have calculated the VDOS of the sodium atoms and the 
atoms of the matrix, by Fourier transforming the corre- 
sponding velocity autocorrelation functions. It is known 
that such a method reproduces only approximately the 
VDOS, since at finite temperature (here T ~ 1900 K), 
the harmonic hypothesis does not fully hold but this ap- 
proximation will not be crucial for our arguments. In 
Fig. 4 (a) the VDOS of the oxygen atoms for diff'erent 
values of fi are represented (a similar picture could be 
drawn for the Si atoms). The VDOS for /x = 1 is of 
course close to the one of amorphous silica obtained ex- 
perimentally except that it is less structured, which is 
a known drawback of the BKS potential [20]. Such a 
VDOS is typical of a broad band of vibrational states, 
extending from = up to i^i ~ 35 THz, as a result 
of the coupling between neighboring oxygen and silicon 
atoms forming a strong random covalent network. When 
increasing the mass of the atoms of the matrix, one ob- 
serves a shrinkage of the VDOS towards low frequencies. 
The top of the oxygen band, f°ax> been plotted ver- 
sus /i in Fig 4 (c): as expected from standard solid state 
physics i^max varies like The situation is quite 

different for the Na VDOS represented in Fig. 4 (b) for 
three typical values of fi: it looks like a well defined peak 
centered at a frequency vf^ close to 5 THz for small 
values of /x. Such a peak can be interpreted as being 
mostly due to the vibrations of the sodium atoms in their 
cage. Using the approximate formula for such a motion, 
R-^y^kT/mj<s^, where i? = 1 A is a typical size of the 
cage deduced from the MSD and rriNa is the mass of a 
sodium atom, one finds a frequency of about 8 THz in 
reasonable agreement with i/f'^. The broadening of the 
peak is then certainly due to the polydispersity of the 
R values as well as to the coupling of the sodium atoms 
with other species. With increasing /x a second peak at 
1/2^ grows out of the high frequency part and increases 
while the peak at uf'^ decreases and shifts towards low 
frequencies. The variation of v^'^ and 1^2^ with n is rep- 
resented in Fig. 4 (c) where the main peak of the Na 
VDOS is represented by the filled symbols. One sees 
that up to /i = 30 the peak at vf^ is the main peak 
while for higher values of /x the peak at fl*^ becomes the 
main peak. Once the transition is fulfilled, the secondary 
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Na peak follows the variation of f^^x with ^ while the 
principal peak remains at a constant frequency. Such a 
behavior is typical of hybridization effects as commonly 
seen in electronic and vibrational systems. It is due to 
the coupling between the sodium atoms and the atoms 
of the matrix. The characteristic value no for which the 
principal peak in the Na VDOS starts to escape from the 
broad oxygen band can be estimated by equating the top 

— 1/2 

of the band ui/jIq and the typical sodium frequency, 
5 THz. This gives /io ~ 50 (indicated by the arrow in 
Fig. 4 (c)), which is remarkably close to the value of 30 
above which we have observed the change in the diffusive 
properties of the Na atoms. 

Therefore the above analysis of the vibrational den- 
sity of states provides a very simple explanation of the 
predominance of channel diffusion for /i smaller than /io- 
In that case the sodium frequency peak lies within the 
limits of the vibrational band of the matrix as shown in 
Fig. 4 (c). Therefore there always exists a vibrational 
mode of the matrix with the same frequency as the one 
of the sodium atoms. Using a simple picture in direct 
space, the vibrational amplitudes of the sodium atoms 
can become very large due to "resonance conditions", 
giving them the possibility to escape from their cage and 
to jump to a neighboring cage, as it is expected in a 
local picture of the diffusive motion. Moreover, at the 
location where a sodium atom is connected to an oxy- 
gen atom via a covalent-like bond, it is clear that the 
local hybridization is important and this explains why 
the location of the channels, i.e. the preferential path- 
ways for the sodium diffusive motion, is strongly corre- 
lated to the position of the non-bridging oxygen atoms. 
On the contrary, for /i > /io, the vibrational peak of the 
sodium atoms is located outside the VDOS of the matrix 
(Fig. 4 (c)) and the resonance conditions are more diffi- 
cult to be fulfilled. As a consequence, the sodium atoms 
stay much longer confined in their cages. This picture is 
consistent with our interpretation of the increasing depth 
of the first minimum of the Na velocity autocorrelation 
function, as being due to an increasing stiffness of the 
internal cage boundaries when is increased above /xq- 

In conclusion, we have demonstrated the essential role 
played by the vibrations of the atoms of the matrix in 
the existence of the channel diffusion of sodium inside 
an amorphous silica matrix. In a sense this is close to 
the concept of "matrix-mediated-coupling" recently used 
to interpret the mixed cation effect in glasses [21] except 
that we show here that the direct hybridization between 
the ionic modes and the modes of the matrix is necessary 
to insure a fast diffusion process of the ions. This result 
is also coherent with previous studies showing indirectly 
that the sodium dynamics should be intimately linked 
to the one of the matrix [15]. Of course it would be 
extremely interesting to extend the present study to other 
kinds of ions, such as Li, H, etc., to test the generality 



of our interpretation (according to preliminary results it 
appears indeed that a similar behavior is observed for Li 
in Si02 [22]). This would constitute a real improvement 
in the understanding of the mechanisms of ionic transport 
in random media. 

We thank Pr. K. Funke and Pr. A. Hcuer for in- 
teresting discussions. Part of the numerical calcula- 
tions were done at the "Centre Informatique National 
de I'Enseignement Superieur" (CINES) in Montpelher. 
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FIG. 1. Plot of the mean square displacement R^{t) of the 
sodium atoms for different values of /i. Insot: Plot of the 
characteristic diffusion time tmsd (see text for definition) as 
a function of //. 
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FIG. 2. Plot of ^, the quantity used to define the channels 
(see text for further details) as a function of /u. 
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FIG. 3. Plot of the velocity autocorrelation function i?(t) 
of the sodium atoms for different values of /j. Inset: Plot of 
the intensity of the first minimum '&rn as a function of //. 
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FIG. 4. (a): Plot of the VDOS of the oxygen atoms for different values of jj,. (b): Plot of the VDOS of the sodium atoms 
respectively for /i = 1, 10, and 30. (c): Plot of the highest frequency of the oxygen VDOS (•), the frequency of the low energy 
peak (v) and the high energy peak (A) of the sodium VDOS as a function of fj,. When these symbols are filled they indicate 
the main peak of the Na VDOS. 
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